Introduccion modelos jerarquicos (I)
1 Introducción
Los modelos jerárquicos (mixtos o multinivel), corresponden a un grupo de modelos estadísticos que permiten incorporar explícitamente la jerarquía de la variación en los datos. Esto es, estiman la distribución posterior de los parámetros del modelo a través agrupamientos parciales (partial pooling) de los grupos, donde los parámetros particulares de cada grupo se estima de una distribución de probabilidad de la “poblacional” de grupos — i.e. una distribución de probabilidad de coeficientes. Dado que todos los grupos provienen de una misma distribución de probabilidad, las observaciones de un grupo ayudan a estimar los coeficientes de otro. Lo que resulta ventajoso cuando el esfuerzo de muestreo no está balanceado entre grupos. Veamos un ejemplo, aclaremos la idea de agrupamientos en los datos y sus implicaciones para la estimación de los coeficientes en un modelo jerárquico.
Modelemos la probabilidad de supervivencia de huevos en 50 nidadas de una ave cualquiera en un parche de bosque. Tenemos tres alternativas:
Modelar los datos asumiendo que todos los nidos tienen igual probabilidad \(p\) de ser depredados (pooling), y por lo tanto estimando un único parámetro común para todos. Esto es, asumimos que la variación entre nidos es cero.
Estimar de manera independiente el coeficiente \(p\) para cada nido (no pooling).
Estimar \(p\) para cada nido, pero asumiendo que cada \(p_i\) hace parte de una distribución de probabilidad \(Binomial(p)\) de la población de nidos (partial pooling).
Simulemos los datos
Ajustemos el modelo pooling.
Modelo no pooling
Modelo partial pooling
par_pool <-
ulam(
alist(
surv ~ dbinom(huevos, p),
logit(p) <- alpha[nido],
alpha[nido] ~ dnorm(mu, sigma), # previa de la población de nidos (previa adaptable) (hiper parámetros)
mu ~ dnorm(5, 2.5), # previa del nido promedio (hiper previa)
sigma ~ dexp(1) # (hiper previa)
), data = dat, chains = 3, cores = 3, iter = 5e3, warmup = 1e3,
log_lik = T
)El modelo partial pooling recibe el nombre de modelo jerárquico, porque la estimación de los parámetros se realiza en varios niveles. Veamos en notación matemática el modelo anterior:
\[ \begin{align} &Sup.~huevos_i \sim Binomial(huevos_i,~p_i) & ,~Funcion~de~likelihhod \\ &logit(p_i) \sim \alpha[nido_i] &,~Modelo~lineal \\ &\alpha[nido] \sim Normal(\mu_j, \sigma_j) &,~Previa~para~la~población~de~nidos \\ &\mu_j \sim Normal(5, 2.5) &,~Previa~para~cada~nido~(hiperparámetro) \\ &\sigma_j \sim exponential(1) &,~Previa~para~cada~nido~(hiperparámetro) \\ \end{align} \] \(\mu_j\) y \(\sigma_j\) se conocen como hiperparámetros porque son parámetros de parámetros (la previa poblacional), por la misma razón \(Normal(5, 2.5)\) y \(exponential(1)\) se les llama hiper-previas.
Usemos el WAIC para comparar los tres modelos en términos de su capacidad predictiva fuera de la muestra:
WAIC SE dWAIC dSE pWAIC weight
par_pool 200.8666 6.298094 0.00000 NA 20.644381 9.999999e-01
pooling 232.7893 17.290137 31.92270 11.732154 2.131153 1.169699e-07
no_pooling 238.3174 5.153536 37.45076 9.454381 38.346678 7.373452e-09El WAI del modelo jerárquico es notablemente más bajo comparado con los modelos pooling y no pooling (i.e. tiene un mayor poder predictivo). Veamos qué tan precisa es la estimación del log-odds de los modelos jerárquico y pooling:
log_odds_P <- extract.samples(no_pooling)
log_odds_parP <- extract.samples(par_pool)
plot(NULL, xlim = c(-5, 5), ylim = c(0, 1.2), xlab = 'log-odds nidos',
ylab = 'Densidad')
for (i in 1:500) {
curve(dnorm(x, log_odds_parP$mu[i], log_odds_parP$sigma[i]),
lwd = 0.1, add = T, col = "red")
lines(density(log_odds_P$alpha[i, ]), lwd = 0.1, col = 'tan1')
}
curve(dnorm(x, mean(log_odds_parP$mu), mean(log_odds_parP$sigma)), add = T,
col = 'red', lwd = 2.5)
lines(density(apply(log_odds_P$alpha, 2, mean)), col = 'tan1', lwd = 2.5)
lines(density(logit(rbeta(1e3, 3, 5))), col = 'black', lwd = 2.5)
Ambos modelos predicen bien el promedio del log-odd real (negro). Sin embargo, notemos la amplitud de las distribuciones posteriores, la estimación del modelo jerárquico (rojo) es mucho más precisa comparada con la del modelo pooling (naranja).
Ahora comparemos el error de la probabilidad de supervivencia estimada por ambos modelos, con la probabilidad “real” — la simulada a partir de \(Beta(\phi_1 = 3, \phi_2= 5)\).
post_parPool <- extract.samples(par_pool)
post_parPool <- as_tibble(post_parPool$alpha)
colnames(post_parPool) <- paste('nido', 1:n, sep = '')
post_parPool <- do.call('rbind', apply(post_parPool, 2, simplify = 'list', FUN =
function(x) {
tibble(li = quantile(x, 0.025),
ls = quantile(x, 0.975),
mu = mean(x),
mod = 'Par pooling')
}))
post_parPool$x <- paste('nido', 1:n, sep = '')
post_noPool <- extract.samples(no_pooling)
post_noPool <- as_tibble(post_noPool$alpha)
colnames(post_noPool) <- paste('nido', 1:n, sep = '')
post_noPool <- do.call('rbind', apply(post_noPool, 2, simplify = 'list', FUN =
function(x) {
tibble(li = quantile(x, 0.025),
ls = quantile(x, 0.975),
mu = mean(x),
mod = 'No pooling')
}))
post_noPool$x <- paste('nido', 1:n, sep = '')
mods <- rbind(post_noPool, post_parPool)
for (i in 1:3) mods[[i]] <- inv_logit(mods[[i]])ggplot() +
geom_errorbar(data = mods, aes(x = x, ymin = li, ymax = ls, color = mod),
position = position_dodge(width = 0.7), width = 0) +
geom_point(data = mods, aes(x = x, y = mu, color = mod),
position = position_dodge(width = 0.7)) +
geom_point(data = tibble(y = p,
x = paste('nido', 1:n, sep = '')),
aes(x, y)) +
geom_hline(yintercept = mean(p), linetype = 2) +
labs(y = 'Probabilidad', x = NULL) +
coord_flip() +
theme_bw() +
theme(panel.grid = element_blank(),
legend.title = element_blank(),
legend.position = 'top')
Los puntos negros corresponden a la probabilidad de supervivencia real. Pareciera que los valores predichos por el modelo jerárquico están, en general, más cercanos a la probabilidad de supervivencia real de los nidos. Calculemos este error de estimación y veamos cómo se comporta con el tamaño de la muestra (i.e. número de huevos observados por nido).
mods$real_p <- rep(p, 2)
mods$sample_size <- c(huevos, huevos)
mods$error <- mods %$% abs(mu - real_p)
ggplot(mods, aes(sample_size, error, color = mod)) +
geom_point() +
geom_smooth(method = 'lm', se = F) +
labs(x = 'Tamaño de la muestra', y = 'Error P(supervivencia)') +
theme_bw() +
theme(panel.grid = element_blank(),
legend.title = element_blank(),
legend.position = 'top')
En ambos modelos el error decrece con el aumento del tamaño de muestra; es natural, más datos por lo general proveen más información y, por lo tanto, una estimación más precisa de los coeficientes. Notemos, sin embargo, que en promedio el modelo jerárquico (par pooling) estima la probabilidad de supervivencia con un menor error, aunque para ambos modelos este valor tiende a converger con el incremento del tamaño de muestra.
2 Modelos con más de un agrupamiento
En el ejemplo anterior modelamos la probabilidad de supervivencia de huevos en 50 nidadas del ave X ubicadas en un parche de bosque cualquiera. Supongamos que la investigadora decidió expandir la investigación, duplicar el número de nidadas analizadas y distribuir su esfuerzo de muestreo en 9 parches de bosques adicionales. Ahora, dado el nuevo diseño, la supervivencia de los huevos no solo varía en función del nido, sino también dado el parche de bosque en el que se encuentra el nido.
Generemos los datos y describamos el modelo en notación matemática:
n_nidos <- 100
n_parches <- 10
set.seed(555)
p <- rbeta(n_nidos, 3, 5)# probabilidad de visita (5 min) sin efecto
set.seed(555)
huevos <- rpois(n_nidos, 10)
set.seed(555)
var_parche <- rnorm(1e3, 0, 0.05)
set.seed(555)
var_parche <- sample(var_parche[var_parche >= 0], n_parches, replace = F)
var_parche <- sapply(var_parche, simplify = 'array', FUN =
function(x) rnorm(n_parches, 0, x))
var_parche <- as.vector(var_parche)
set.seed(555)
huevos_obs <- rbinom(n_nidos,
prob = p + var_parche,
size = huevos)
dat <-
list(
n = huevos,
superv = huevos_obs,
parche = rep(1:10, each = 10),
nido = 1:100
)\
\[ \begin{align} &Sup.~huevos_i \sim Binomial(huevos_i,~p_i) & ,~Funcion~de~likelihhod \\ &logit(p_i) \sim \alpha[nido_i] + \psi[parche] &,~Modelo~lineal \\ &\psi[parche] \sim Normal(0, \tau) &,~Previa variación~entre~parches\\ &\tau \sim exponential(1) &, ~ Hiper-parámetro/previa~de~la~previa~del~parche\\ &\alpha[nido] \sim Normal(\mu_j, \sigma_j) &,~Previa~para~la~población~de~nidos \\ &\mu_j \sim Normal(5, 2.5) &,~Previa~para~cada~nido~(hiperparámetro) \\ &\sigma_j \sim exponential(1) &,~Previa~para~cada~nido~(hiperparámetro) \\ \end{align} \]
El parche corresponde a un factor que altera la supervivencia de los huevos. Es decir, \(p_i\) reduce o aumenta en función del parche. Dichas oscilaciones las incorporamos en el modelo en forma de una distribución \(Normal(0, \tau)\), donde \(\tau\) controla la ‘magnitud’ de los cambios positivos o negativos.
Ajustemos el modelo:
mod <-
ulam(
alist(
superv ~ dbinom(n, p),
logit(p) <- alpha[nido] + psi[parche],
# previas poblacionales
alpha[nido] ~ dnorm(mu_j, alpha_j),
psi[parche] ~ dnorm(0, tau),
# Hiperparámetros (previas nido y parche)
mu_j ~ dnorm(5, 2.5),
alpha_j ~ dexp(1),
tau ~ dexp(1)
), data = dat, chains = 3, warmup = 500, iter = 6e3, cores = 3
) mean sd 5.5% 94.5% n_eff Rhat4
mu_j -0.6909097 0.1252373 -0.89358515 -0.4977533 7983.075 0.9999535
alpha_j 0.8148713 0.1180152 0.63479064 1.0067738 3863.767 0.9999266
tau 0.1327899 0.1041278 0.02055271 0.3232691 1365.567 1.0029164Tenemos un parámetro por cada nido, uno por cada parche y su parámetro de dispersión y la probabilidad de supervivencia promedio. Antes de discutir las implicaciones de los parámetros, grafiquemos las distribuciones predictivas posteriores y estimemos la precisión del modelo para representar las observaciones.
post <- extract.samples(mod)
post <- lapply(post, FUN =
function(x) {
apply(as.matrix(x), 2, inv_logit)
})
ppcheck <- sapply(1:100, simplify = 'array', FUN =
function(x) rbinom(1e3,
size = dat$n,
prob = post$mu_j[i]))
plot(NULL, xlim = c(0, 10), ylim = c(0, 0.6))
for (i in 1:100) lines(density(ppcheck[i, ]), lwd = 0.15, col = 'blue')
lines(density(dat$superv), col = 'red', lwd = 3)
abline(v = c(mean(apply(ppcheck, 2, mean)), mean(dat$superv)),
col = c('blue', 'red'), lty = 3)
El modelo se ajusta adecuadamente a los datos, incluso el promedio estimado de huevos supervivientes es básicamente igual al observado.
Veamos las implicaciones de los parámetros
par(mfrow = c(1, 2))
plot(NULL, xlim = c(0, 1), ylim = c(0, 6), xlab = expression(italic('P')['nido']),
ylab = 'Densidad', main = 'Nido')
for (i in 1:100) lines(density(post$alpha[, i, drop = T]), lwd = 0.2)
lines(density(apply(post$alpha, 2, mean)), lwd = 3, col = 'red')
lines(density(p + var_parche), lwd = 3, col = 'blue')
plot(density(post$alpha_j), main = "Parche", xlim = c(0.47, 0.81),
xlab = expression(sigma[italic(P)]), ylab = '')
lines(density(post$tau), lty = 2)
Izquierda. Cada línea de densidad negra corresponde a la distribución posterior de \(p\) para cada nido, la línea roja denota el promedio entre nidos —igual que el parámetro \(\mu_i\) en el modelo, y la línea azul corresponde a la probabilidad real (i.e. probabilidad de supervivencia, \(Beta(3, 5)\), más variación dado el parche, \(Normal(0, \tau)\)). Derecha. Distribución posterior de los parámetros de dispersión de los nidos (línea sólida) y los parches (línea discontinua).
La probabilidad de supervivencia promedio es del ~40%.
De hecho, el modelo fue capaz de recuperar la distribución de probabilidad generadora de los datos (línea azul). También, con seguridad, podemos afirmar que la variación en la probabilidad de supervivencia entre nidos es mayor a la producida por los diferentes parches. Podríamos generar un nuevo modelo sin el factor parche, compararlos con WAIC, y evaluar el aporte del parche en la capacidad predictiva del modelo.
Podríamos concluir (ponele…), que los factores que determinan la probabilidad de supervivencia de huevos en nidos del ave X operan principalmente a escala local —la principal variación proviene del nido, mientras que el aporte del parche es tangencial. Supongamos que la investigadora alcanzó la misma conclusión y, en un intento de identificar dichos factores locales moduladores de la supervivencia, decidió medir la proporción de cobertura vegetal en un radio de 1 m alrededor de cada nido. Veamos:
3 Transisiones divergentes y reparametrización
De nuevo simulemos los datos y definamos el modelo.
n_nidos <- 100
n_parches <- 10
set.seed(555)
p <- rbeta(n_nidos, 3, 5)
set.seed(555)
huevos <- rpois(n_nidos, 10)
set.seed(555)
var_parche <- rnorm(1e3, 0, 0.1)
set.seed(555)
var_parche <- sample(var_parche[var_parche >= 0], n_parches, replace = F)
var_parche <- sapply(var_parche, simplify = 'array', FUN =
function(x) rnorm(n_parches, 0, x))
var_parche <- as.vector(var_parche)
tot_p <- p + var_parche
set.seed(555)
pro_veg <- rbeta(1e3, 2, 1.5)
temp2 <- quantile(pro_veg, probs = seq(0, 1, by = 0.2))
temp <- quantile(tot_p, probs = seq(0, 1, by = 0.2))
veg <- vector('double', 100)
set.seed(555)
for (i in 1:(length(temp)-1)) {
veg[which(tot_p >= temp[i] & tot_p <= temp[i+1])] <-
sample(pro_veg[pro_veg > temp2[i] & pro_veg <= temp2[i+1]],
size = length(which(tot_p >= temp[i] & tot_p <= temp[i+1])),
replace = F) +
rnorm(length(which(tot_p >= temp[i] & tot_p <= temp[i+1])),
0, 0.05)
}
set.seed(555)
beta_veg <- 0.1
set.seed(555)
huevos_obs <- rbinom(n_nidos,
prob = p + var_parche + beta_veg*veg,
size = huevos)
dat <-
list(
n = huevos,
superv = huevos_obs,
parche = rep(1:10, each = 10),
prop_veg = veg,
nido = 1:100
)En principio, ajustemos el modelo anterior, pero incluyamos un parámetro \(\beta\) que controle la relación entre la supervivencia de huevos y la proporción de cobertura vegetal. Esto es:
\[ \begin{align} &Supervivencia_i \sim Binomial(huevos_i,~p_i) \\ &logit(p_i) \sim \alpha[nido_i] + \psi[parche] + \beta\times{prop. vegetación}\\ &\beta \sim Normal(0.5, 1)\\ &\psi[parche] \sim Normal(0, \tau) \\ &\tau \sim exponential(1) \\ &\alpha[nido] \sim Normal(\mu_j, \sigma_j) \\ &\mu_j \sim Normal(5, 2.5) \\ &\sigma_j \sim exponential(1) \\ \end{align} \]
mod2.1 <-
ulam(
alist(
superv ~ dbinom(n, p),
logit(p) <- alpha[nido] + psi[parche] + beta*prop_veg,
beta ~ dnorm(0.5, 1),
alpha[nido] ~ dnorm(mu, sigma),
psi[parche] ~ dnorm(0, tau),
mu ~ dnorm(10, 5),
sigma ~ dexp(1),
tau ~ dexp(1)
), data = dat, iter = 5e3, chains = 3, cores = 3, warmup = 500, log_lik = T)Warning: There were 167 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
https://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemsWarning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-essWarning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess mean sd 5.5% 94.5% n_eff Rhat4
beta 3.595314476 0.35928991 3.01816575 4.1765160 526.65200 1.0010461
alpha[1] -2.457878817 0.40454544 -3.12583061 -1.8188119 1339.42627 1.0002388
alpha[2] -2.277407787 0.40944078 -2.91516253 -1.5970097 1210.61577 1.0005329
alpha[3] -2.508981041 0.40497124 -3.16766775 -1.8812802 1173.58039 1.0013922
alpha[4] -1.794322959 0.48447986 -2.49228771 -0.9665070 356.62648 1.0078461
alpha[5] -2.519908983 0.44351845 -3.25540175 -1.8268143 1270.81367 1.0011101
alpha[6] -2.577136666 0.43061569 -3.28364708 -1.9262652 1300.51055 1.0013036
alpha[7] -2.760709558 0.42216926 -3.46121789 -2.1254098 530.32271 1.0044241
alpha[8] -2.467119598 0.41499869 -3.14923220 -1.8289025 1442.23690 1.0009015
alpha[9] -2.704083490 0.42313414 -3.41905010 -2.0801254 748.73570 1.0036718
alpha[10] -2.301928092 0.40441078 -2.94755552 -1.6484585 1379.54916 1.0000480
alpha[11] -2.229128931 0.41844135 -2.88242800 -1.5387038 1440.41597 1.0012571
alpha[12] -2.411363697 0.41444796 -3.08411751 -1.7607217 1553.08090 1.0003774
alpha[13] -2.607945297 0.41058722 -3.29292766 -1.9948955 1055.70379 1.0025925
alpha[14] -2.431540108 0.41283503 -3.10200291 -1.7764155 1511.11343 1.0002736
alpha[15] -2.052812280 0.42064947 -2.68342129 -1.3434536 783.40764 1.0026807
alpha[16] -2.493587270 0.41302485 -3.16980443 -1.8560529 1168.97450 1.0009772
alpha[17] -2.458187564 0.42133131 -3.14801694 -1.8124054 1450.59861 1.0006936
alpha[18] -2.621045350 0.40953265 -3.29505827 -1.9995087 843.60234 1.0022919
alpha[19] -2.353543362 0.40278911 -2.99479003 -1.7101148 1237.95999 1.0005130
alpha[20] -2.657643595 0.44549940 -3.41497858 -1.9909879 1006.20579 1.0022819
alpha[21] -2.321527448 0.38667782 -2.92962614 -1.6946909 1150.41639 1.0005141
alpha[22] -2.629109697 0.43343928 -3.36354857 -1.9837147 1129.13437 1.0016718
alpha[23] -2.180726014 0.43734074 -2.84959102 -1.4500614 1218.82037 1.0012843
alpha[24] -2.614490628 0.39559592 -3.27618247 -2.0181911 958.43485 1.0011360
alpha[25] -2.650393188 0.42438513 -3.36241883 -2.0192037 1008.77310 1.0020629
alpha[26] -2.365461248 0.40718732 -3.01900304 -1.7142745 1321.46882 1.0006008
alpha[27] -2.509948626 0.39946661 -3.16626286 -1.8967537 1239.43138 1.0003394
alpha[28] -2.278214844 0.40093541 -2.90522426 -1.6193673 1329.27683 1.0003773
alpha[29] -2.567865861 0.41369875 -3.24457902 -1.9282573 1160.19479 1.0008173
alpha[30] -2.478632276 0.40304652 -3.13810077 -1.8544772 1500.89373 1.0007618
alpha[31] -2.312766734 0.41047348 -2.95501201 -1.6449861 1506.85865 1.0005722
alpha[32] -2.569041524 0.40656225 -3.25271831 -1.9513771 1049.58268 1.0019372
alpha[33] -2.326444078 0.41076908 -2.97731856 -1.6579664 1412.57380 1.0004592
alpha[34] -2.469530743 0.41862913 -3.15113707 -1.8224302 1418.88484 1.0008294
alpha[35] -2.548925341 0.38684052 -3.17747621 -1.9388834 1046.26018 1.0017671
alpha[36] -2.247293524 0.41770651 -2.90318151 -1.5626794 1161.15988 1.0005195
alpha[37] -2.597028490 0.40335690 -3.28605849 -1.9870746 1077.64663 1.0015516
alpha[38] -2.563606912 0.41798973 -3.26475146 -1.9253912 1251.68847 1.0008688
alpha[39] -2.483819210 0.40531024 -3.14361800 -1.8509088 1159.05810 1.0009345
alpha[40] -2.450482041 0.40834690 -3.10967554 -1.8105691 1384.28978 1.0005192
alpha[41] -2.436679092 0.37866398 -3.05387488 -1.8344145 1270.82811 1.0006126
alpha[42] -2.428891132 0.38696800 -3.04878918 -1.8081557 1234.60852 1.0004934
alpha[43] -2.484094143 0.42196198 -3.17000391 -1.8305115 1222.96583 1.0006563
alpha[44] -2.314723400 0.40822952 -2.95748094 -1.6596120 1546.08250 1.0005958
alpha[45] -2.328149559 0.41166604 -2.97816929 -1.6607076 1191.38074 1.0008907
alpha[46] -2.110388976 0.42747284 -2.74667630 -1.3953051 885.77483 1.0019788
alpha[47] -2.220417043 0.42085909 -2.86576509 -1.5237134 1235.02675 1.0013818
alpha[48] -2.575409779 0.41086657 -3.25246285 -1.9491967 1099.48756 1.0016146
alpha[49] -2.580024437 0.44554758 -3.31862080 -1.9121845 1372.52092 1.0016783
alpha[50] -2.252869849 0.44060375 -2.93406829 -1.5322163 1250.24497 1.0010342
alpha[51] -2.329468867 0.38047531 -2.93415676 -1.7195508 1426.60188 1.0000758
alpha[52] -2.552228490 0.40452887 -3.23860640 -1.9388931 1346.90624 1.0011701
alpha[53] -2.302937268 0.42605258 -2.97235465 -1.6139561 1433.33818 1.0005563
alpha[54] -2.193732560 0.40327509 -2.81236385 -1.5269946 1127.71027 1.0009495
alpha[55] -2.428720170 0.43004748 -3.13029678 -1.7539294 1291.63317 1.0006993
alpha[56] -2.464446177 0.42094331 -3.15490255 -1.8178070 1538.93083 1.0007766
alpha[57] -2.364734557 0.44135861 -3.07218988 -1.6585156 1380.35772 1.0003189
alpha[58] -2.242037959 0.43579228 -2.90939803 -1.5167838 1193.37766 1.0010606
alpha[59] -2.249913022 0.39254630 -2.86931747 -1.6184604 1192.21626 1.0005833
alpha[60] -2.340869147 0.38131402 -2.94995276 -1.7185856 1244.84260 1.0004960
alpha[61] -2.402925756 0.42232728 -3.08274580 -1.7269603 1652.15290 1.0001200
alpha[62] -2.700513009 0.42971601 -3.43655789 -2.0758143 922.28764 1.0029325
alpha[63] -2.433434960 0.39821484 -3.07872487 -1.7958529 1434.10741 1.0008657
alpha[64] -2.323375432 0.41568900 -2.97442515 -1.6511621 1229.25043 1.0007165
alpha[65] -2.394897220 0.38693395 -3.01978823 -1.7734423 1312.87883 1.0002364
alpha[66] -2.330150405 0.39374023 -2.96197738 -1.6901639 1660.87544 1.0006027
alpha[67] -2.460155467 0.39599539 -3.10249112 -1.8372686 1284.66227 1.0004768
alpha[68] -2.372065251 0.38544123 -2.98520710 -1.7597650 1212.23438 1.0003598
alpha[69] -2.420746829 0.40519100 -3.07714004 -1.7860632 1430.60209 1.0004190
alpha[70] -2.371236767 0.43149842 -3.05033305 -1.6794719 1332.09496 1.0004548
alpha[71] -2.420274231 0.42997911 -3.12000613 -1.7390683 1496.68057 1.0002610
alpha[72] -2.691577178 0.43810182 -3.43582263 -2.0533363 976.23197 1.0025656
alpha[73] -2.306767830 0.44009674 -2.99502749 -1.5811897 1344.87602 1.0003833
alpha[74] -2.459771936 0.39810709 -3.11135883 -1.8363717 1162.73324 1.0004946
alpha[75] -2.346442994 0.40310483 -3.00056537 -1.7138491 1712.18395 1.0000243
alpha[76] -2.470705116 0.40066453 -3.12399821 -1.8446621 1166.84382 1.0002918
alpha[77] -2.280246378 0.43026652 -2.94779315 -1.5708748 1490.13696 1.0004511
alpha[78] -2.402671772 0.40569179 -3.05762492 -1.7612763 1649.04400 1.0006432
alpha[79] -2.559977221 0.40305572 -3.23906786 -1.9484786 1165.87121 1.0014279
alpha[80] -2.323143103 0.42597150 -2.99714445 -1.6266592 1563.26293 1.0005279
alpha[81] -1.972923079 0.43866971 -2.61405413 -1.2332942 639.78323 1.0038332
alpha[82] -2.110485790 0.41282054 -2.74106448 -1.4234810 1223.65810 1.0016532
alpha[83] -2.777946103 0.45112289 -3.52649238 -2.1086022 543.01421 1.0051363
alpha[84] -2.221864879 0.39547806 -2.83254161 -1.5734623 1156.83471 1.0013480
alpha[85] -2.592925523 0.41009754 -3.27476286 -1.9689319 993.34504 1.0017750
alpha[86] -2.573364878 0.41659644 -3.26273381 -1.9409763 1021.95870 1.0015207
alpha[87] -2.504144856 0.43381565 -3.21785922 -1.8405229 1419.17390 1.0013180
alpha[88] -2.346434955 0.41988744 -3.01199268 -1.6613746 1300.54750 1.0005423
alpha[89] -2.672342285 0.44872690 -3.41957753 -2.0124117 882.50764 1.0032987
alpha[90] -2.356026448 0.38271820 -2.97401388 -1.7401103 1606.41599 1.0001527
alpha[91] -2.477691948 0.43029635 -3.17588727 -1.8098188 1187.56200 1.0006907
alpha[92] -2.406044494 0.41710039 -3.07610502 -1.7431004 1318.51428 1.0001924
alpha[93] -2.520102278 0.41593733 -3.20806799 -1.8792811 1392.41001 1.0006421
alpha[94] -2.700048205 0.42915577 -3.42341215 -2.0590099 654.51129 1.0034368
alpha[95] -2.491952540 0.41342299 -3.15874856 -1.8616216 1377.23927 1.0006058
alpha[96] -2.253950269 0.39357585 -2.86982542 -1.6167613 1449.99502 1.0006017
alpha[97] -1.900434166 0.43153642 -2.53358118 -1.1731612 461.10270 1.0058978
alpha[98] -2.308064967 0.41549455 -2.95960922 -1.6286538 1235.82223 1.0006481
alpha[99] -2.017552109 0.46549170 -2.69204676 -1.2177935 783.29631 1.0024455
alpha[100] -2.527495401 0.41351430 -3.20901590 -1.8877412 1185.05478 1.0014096
psi[1] -0.025775678 0.11544014 -0.22432115 0.1392797 3294.25620 0.9999892
psi[2] -0.018552110 0.11730276 -0.21568062 0.1505488 2785.25535 1.0017538
psi[3] -0.049780327 0.12320614 -0.27491257 0.1058420 2057.53113 1.0008003
psi[4] -0.042853466 0.12143647 -0.26470359 0.1130292 2580.78784 1.0006390
psi[5] 0.037560248 0.12368562 -0.13265734 0.2498864 1505.66122 1.0028878
psi[6] 0.071820506 0.13903003 -0.08477647 0.3431003 772.00399 1.0024765
psi[7] -0.007872746 0.11453627 -0.19301812 0.1711340 4150.65377 1.0006695
psi[8] -0.010399873 0.12001716 -0.20276708 0.1671643 2575.57454 1.0001531
psi[9] -0.003266585 0.11271163 -0.18143152 0.1729741 4887.76039 1.0013559
psi[10] 0.050448667 0.12500917 -0.10248053 0.2825607 2490.42487 1.0010065
mu -2.411953712 0.22701454 -2.79049528 -2.0527624 448.91918 1.0013587
sigma 0.378304533 0.13104980 0.14881660 0.5796904 99.04207 1.0266395
tau 0.117156955 0.09331631 0.01332070 0.2862240 548.30054 1.0063168[1] 167El modelo tiene un un número tal de transiciones divergente que el effective sampling size de algunos parámetros es inferior a 1000 y por lo tanto no deberíamos estar muy confiados de esta estimación. Esto es común en modelos jerárquicos, y sucede porque el algoritmo explora incorrectamente algunas regiones del espacio de los parámetros. Podríamos intentar solucionarlo de varias formas: (i) incrementando el número de iteraciones (no es lo ideal porque no lidiaríamos con el problema de eficiencia en la exploración del algoritmo); (ii) incrementando la tasa de aceptación de valores que explora el algoritmo controlando los “saltos” entre uno y otro valor estimado (leapfrog step); (iii) reparametrización del modelo. Probemos la alternativa ii incrementando el adapt_delta = 0.99:
mod2.2 <-
ulam(
alist(
superv ~ dbinom(n, p),
logit(p) <- alpha[nido] + psi[parche] + beta*prop_veg,
beta ~ dnorm(0.5, 1),
alpha[nido] ~ dnorm(mu, sigma),
psi[parche] ~ dnorm(0, tau),
mu ~ dnorm(10, 5),
sigma ~ dexp(1),
tau ~ dexp(1)
), data = dat, iter = 5e3, chains = 3, cores = 3, warmup = 500, log_lik = T,
control=list(adapt_delta=0.99))Warning: There were 9 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: There were 4107 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceededWarning: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
https://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemsWarning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-essWarning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess mean sd 5.5% 94.5% n_eff Rhat4
beta 3.534573747 0.37259298 2.95886020 4.15473123 328.4064 1.005028
alpha[1] -2.416412118 0.41051705 -3.08806560 -1.77532937 982.4852 1.001479
alpha[2] -2.229859343 0.42373906 -2.89582640 -1.54397681 780.6298 1.003327
alpha[3] -2.467617597 0.41117334 -3.13988618 -1.82562336 825.0201 1.000365
alpha[4] -1.740585592 0.48536503 -2.47283239 -0.92985372 422.7919 1.012870
alpha[5] -2.476273993 0.44956415 -3.20874566 -1.78523204 961.0587 1.000946
alpha[6] -2.531834668 0.44606735 -3.26590083 -1.86166994 1076.3537 1.000382
alpha[7] -2.716679015 0.42672097 -3.42591248 -2.07359751 714.8002 1.001005
alpha[8] -2.428283539 0.41514553 -3.10786505 -1.79419590 1158.3512 1.000919
alpha[9] -2.665199335 0.42596714 -3.37363883 -2.02730384 817.2138 1.000604
alpha[10] -2.260465045 0.40741377 -2.90624137 -1.60671166 864.0465 1.002952
alpha[11] -2.196319423 0.42893815 -2.87632124 -1.50515823 898.1104 1.003725
alpha[12] -2.373543769 0.42941373 -3.06907828 -1.70539161 1119.8227 1.000706
alpha[13] -2.567940501 0.41101015 -3.23792888 -1.94904928 1033.0211 1.000295
alpha[14] -2.394741393 0.42586720 -3.06799605 -1.72803592 965.1057 1.001845
alpha[15] -2.013000493 0.43334661 -2.67253539 -1.28417963 608.4973 1.007406
alpha[16] -2.444664055 0.42348527 -3.12595423 -1.78441289 856.7383 1.001759
alpha[17] -2.417887302 0.42279155 -3.10677510 -1.75827056 997.7196 1.001273
alpha[18] -2.590645420 0.41339892 -3.27227192 -1.96328635 953.1719 1.000990
alpha[19] -2.309638591 0.40419450 -2.94755313 -1.67233698 815.4909 1.001976
alpha[20] -2.612163072 0.45114047 -3.36774605 -1.94262232 803.2785 1.000565
alpha[21] -2.278175932 0.39713185 -2.91576777 -1.64197681 756.4174 1.003148
alpha[22] -2.592542333 0.43882064 -3.32457686 -1.94077608 987.2774 1.000198
alpha[23] -2.132499658 0.44892962 -2.81197547 -1.38848274 693.4272 1.005049
alpha[24] -2.572044584 0.39808776 -3.22143370 -1.96679445 817.8178 1.000319
alpha[25] -2.608586231 0.42670480 -3.31919293 -1.96332699 823.5024 1.000469
alpha[26] -2.323680456 0.41414945 -2.98529727 -1.67704229 827.2137 1.001435
alpha[27] -2.466300417 0.40302311 -3.12356867 -1.85140215 897.8105 1.000506
alpha[28] -2.239316078 0.41500474 -2.89770144 -1.55923644 771.2270 1.002240
alpha[29] -2.517838255 0.41866550 -3.21189061 -1.88458470 864.8742 1.000530
alpha[30] -2.438864142 0.40607820 -3.10079098 -1.81558668 989.3471 1.000970
alpha[31] -2.269959945 0.41765211 -2.93144948 -1.60539864 894.3477 1.002004
alpha[32] -2.525270148 0.40414121 -3.18213726 -1.90996603 781.7237 1.000765
alpha[33] -2.283052306 0.41505959 -2.94821072 -1.62233997 953.7765 1.002552
alpha[34] -2.426785875 0.42118719 -3.11129915 -1.77065861 922.6412 1.000513
alpha[35] -2.506826121 0.39765635 -3.15941182 -1.89605702 775.2031 1.000723
alpha[36] -2.200848445 0.41867226 -2.85132188 -1.51833509 656.8891 1.004631
alpha[37] -2.560067561 0.40954725 -3.23743742 -1.94904738 870.6989 1.001086
alpha[38] -2.522574093 0.42374036 -3.21851255 -1.88054402 977.1046 1.000601
alpha[39] -2.426387369 0.41521324 -3.09702547 -1.78058852 806.4735 1.001198
alpha[40] -2.410004273 0.41374343 -3.08070633 -1.77362056 1037.7660 1.001060
alpha[41] -2.398258237 0.38342446 -3.02565005 -1.79791648 877.6827 1.001061
alpha[42] -2.379299991 0.39427141 -3.01293826 -1.76091074 846.8422 1.001943
alpha[43] -2.439558280 0.43264369 -3.13793012 -1.76764886 817.1895 1.001411
alpha[44] -2.278727957 0.42267671 -2.95213202 -1.60916671 1115.5322 1.001588
alpha[45] -2.286104018 0.42179039 -2.95830013 -1.60645270 769.4482 1.003044
alpha[46] -2.063540445 0.43659464 -2.73551440 -1.33566491 640.9005 1.006017
alpha[47] -2.173501278 0.41693084 -2.81572463 -1.48979513 777.5116 1.004163
alpha[48] -2.530892649 0.41783035 -3.23032247 -1.89884808 890.3687 1.000199
alpha[49] -2.531495689 0.43995766 -3.26124147 -1.86066801 1068.5959 1.000512
alpha[50] -2.206669652 0.44118857 -2.88892948 -1.47752314 784.4798 1.003320
alpha[51] -2.293453896 0.38631711 -2.91483272 -1.67374945 980.8741 1.002378
alpha[52] -2.513780464 0.41033036 -3.20374208 -1.89966426 1138.8731 1.000729
alpha[53] -2.264403459 0.42980812 -2.95398046 -1.57619537 931.2920 1.002352
alpha[54] -2.146999355 0.41379377 -2.79545273 -1.47381226 707.7770 1.005053
alpha[55] -2.386708148 0.43173312 -3.07068302 -1.71290654 939.6336 1.001059
alpha[56] -2.431466652 0.42184816 -3.11871348 -1.77674621 1148.7935 1.000474
alpha[57] -2.314998776 0.44570371 -3.00933465 -1.60429596 904.0197 1.002648
alpha[58] -2.192406540 0.43536989 -2.86763545 -1.47450798 663.0579 1.003821
alpha[59] -2.208842952 0.40300971 -2.84405746 -1.55489511 767.7687 1.004837
alpha[60] -2.304082513 0.38632226 -2.92086791 -1.69540451 818.9684 1.002934
alpha[61] -2.366794267 0.42920002 -3.05823148 -1.69583780 1044.6958 1.001331
alpha[62] -2.664590338 0.44082439 -3.41376400 -2.03025525 939.9558 1.000524
alpha[63] -2.390556550 0.40398519 -3.04415002 -1.75676191 880.6805 1.001221
alpha[64] -2.278061717 0.41966724 -2.94480884 -1.60227033 724.1856 1.003716
alpha[65] -2.350699916 0.39997057 -2.99940961 -1.73365618 831.9393 1.001947
alpha[66] -2.286650487 0.39419693 -2.91218153 -1.65556391 998.1956 1.001960
alpha[67] -2.418392648 0.41631218 -3.09760189 -1.77276231 996.7079 1.000993
alpha[68] -2.323437292 0.39171900 -2.95470262 -1.70816469 742.7677 1.001585
alpha[69] -2.377909003 0.41321068 -3.03781176 -1.72939794 869.8058 1.001442
alpha[70] -2.320081040 0.43170502 -3.00400121 -1.64496513 857.0707 1.002286
alpha[71] -2.375547989 0.43843739 -3.08202563 -1.69099891 930.7919 1.001226
alpha[72] -2.648385368 0.42980678 -3.37462540 -2.00708652 911.8844 1.001252
alpha[73] -2.261734282 0.44559849 -2.96010342 -1.54003187 859.0674 1.002169
alpha[74] -2.415081922 0.41397685 -3.08288805 -1.76268083 777.6127 1.001164
alpha[75] -2.303437776 0.40815642 -2.96633910 -1.65430499 1096.1020 1.002458
alpha[76] -2.426034456 0.41071982 -3.08492594 -1.79494818 845.8421 1.001021
alpha[77] -2.241938978 0.42743746 -2.90719497 -1.54330824 889.9891 1.003329
alpha[78] -2.366991439 0.41018420 -3.03506813 -1.71650534 1037.6251 1.001170
alpha[79] -2.517145679 0.40815266 -3.19340345 -1.90320143 940.0065 1.000908
alpha[80] -2.270943311 0.43489563 -2.95647413 -1.57328930 836.8444 1.001941
alpha[81] -1.926789021 0.44357364 -2.59718871 -1.18196858 545.1571 1.008878
alpha[82] -2.071993081 0.42073815 -2.72921591 -1.37913712 753.5317 1.006321
alpha[83] -2.733119454 0.44899415 -3.49285913 -2.05622115 726.6571 1.001065
alpha[84] -2.177864810 0.40081826 -2.80656474 -1.51908305 765.3641 1.004211
alpha[85] -2.544169919 0.41762683 -3.23174240 -1.90156265 730.6616 1.000736
alpha[86] -2.527496624 0.42175001 -3.22281404 -1.88351652 861.6106 1.000545
alpha[87] -2.463375460 0.43062165 -3.17576110 -1.80695607 1048.3161 1.000438
alpha[88] -2.300515511 0.42899203 -2.98591063 -1.61457144 827.5832 1.002627
alpha[89] -2.630522526 0.45686011 -3.39865572 -1.94818295 899.4853 1.000471
alpha[90] -2.320892165 0.39121059 -2.94782770 -1.70834621 1067.0629 1.001970
alpha[91] -2.427684067 0.43202844 -3.12728715 -1.75105171 861.3449 1.000977
alpha[92] -2.364591332 0.42663333 -3.04726554 -1.68271675 854.5708 1.001666
alpha[93] -2.483063629 0.42263882 -3.16455828 -1.83529977 1037.1222 1.000556
alpha[94] -2.659450724 0.43147587 -3.37605308 -2.01717652 803.6148 1.001074
alpha[95] -2.449390000 0.41612331 -3.12793739 -1.81293738 998.7187 1.001022
alpha[96] -2.212048552 0.39781034 -2.85109718 -1.57134559 844.5224 1.003144
alpha[97] -1.854978896 0.43459104 -2.52614610 -1.13477149 482.8347 1.010870
alpha[98] -2.264296666 0.41623937 -2.93336925 -1.60061978 779.1146 1.002559
alpha[99] -1.976601938 0.46145522 -2.67580465 -1.19595364 619.3899 1.007477
alpha[100] -2.480719339 0.41669218 -3.15487046 -1.83981472 899.2438 1.000551
psi[1] -0.025234626 0.11368939 -0.22610507 0.13599772 4009.2555 1.000467
psi[2] -0.021863656 0.11748253 -0.22169857 0.14616246 2870.4011 1.000568
psi[3] -0.049378376 0.12183437 -0.27538887 0.09925249 2001.8660 1.000445
psi[4] -0.046631538 0.12316777 -0.26915569 0.10720429 1928.5466 1.000716
psi[5] 0.036370892 0.12003005 -0.12501978 0.25029020 2973.3599 1.000942
psi[6] 0.063622505 0.13185094 -0.08950439 0.30802834 1601.4004 1.001743
psi[7] -0.014648404 0.11620826 -0.21302299 0.15872240 2790.9792 1.000257
psi[8] -0.018573292 0.11809024 -0.21649847 0.15082574 2720.9240 1.000624
psi[9] -0.004411819 0.11296728 -0.18747607 0.17343431 3842.6562 1.000113
psi[10] 0.051254787 0.12613931 -0.10746717 0.28355695 2201.3307 1.000525
mu -2.369202141 0.23826279 -2.77079815 -2.00502114 302.6965 1.004622
sigma 0.381254743 0.12265378 0.18139107 0.57444470 255.0376 1.023096
tau 0.116381663 0.09207426 0.01102023 0.28503742 630.8387 1.004048[1] 9Bien, tenemos menos transiciones divergentes, pero el effective sampling size continua bajo para algunos coeficientes. Probemos entonces reparametrizando el modelo, específicamente una reparametrización no centrada (non-centered reparametrization). Esta consiste en una estandarización que permite que el algoritmo explore el espacio de los parámetros de manera más eficiente. Veamos:
mod2.3 <-
ulam(
alist(
superv ~ dbinom(n, p),
logit(p) <- mu + z[nido] + sigma + # se incluyen en el modelo lineal
x[parche] + tau +
beta*prop_veg,
z[nido] ~ dnorm(0, 1), # intercepto estandarizado del nido
x[parche] ~ dnorm(0, 1), # intercepto estandarizado del parche
beta ~ dnorm(0.5, 1),
mu ~ dnorm(10, 5),
sigma ~ dexp(1),
tau ~ dexp(1),
gq> vector[nido]:alpha <<- mu + z*sigma, # des-estanarización
gq> vector[parche]:psi <<- x*tau # des-estanarización
), data = dat, chains = 3, cores = 3, iter = 5e3, warmup = 500, log_lik = T
)Warning: There were 28 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceededWarning: Examine the pairs() plot to diagnose sampling problems[1] 28 mean sd 5.5% 94.5% n_eff Rhat4
beta 3.39 0.50 2.59 4.20 7147.66 1
mu -3.56 0.99 -5.37 -2.25 4842.66 1
sigma 0.64 0.63 0.04 1.87 8281.73 1
tau 0.64 0.63 0.04 1.87 7679.23 1
psi[1] -0.11 0.46 -0.82 0.42 5609.20 1
psi[2] -0.13 0.46 -0.85 0.36 5268.10 1
psi[3] -0.17 0.46 -0.93 0.31 5265.87 1
psi[4] -0.17 0.45 -0.91 0.30 5783.96 1
psi[5] 0.10 0.44 -0.42 0.78 4427.15 1
psi[6] 0.19 0.49 -0.28 0.95 5152.29 1
psi[7] -0.08 0.43 -0.73 0.45 4826.01 1
psi[8] -0.06 0.44 -0.70 0.49 5721.02 1
psi[9] -0.05 0.43 -0.69 0.50 5818.72 1
psi[10] 0.13 0.44 -0.36 0.85 5138.60 1
alpha[1] -3.64 1.18 -5.78 -2.21 5705.35 1
alpha[2] -3.23 1.01 -5.00 -1.91 5403.13 1
alpha[3] -3.65 1.17 -5.78 -2.21 5236.37 1
alpha[4] -2.32 1.14 -3.98 -0.52 8749.13 1
alpha[5] -3.77 1.32 -6.13 -2.21 5976.82 1
alpha[6] -4.03 1.51 -6.74 -2.32 6041.66 1
alpha[7] -4.03 1.41 -6.65 -2.36 5857.85 1
alpha[8] -3.71 1.25 -5.96 -2.23 6134.18 1
alpha[9] -4.04 1.43 -6.67 -2.35 5766.66 1
alpha[10] -3.29 1.01 -5.12 -1.99 5853.55 1
alpha[11] -3.14 0.99 -4.90 -1.85 6280.25 1
alpha[12] -3.59 1.19 -5.72 -2.16 5999.80 1
alpha[13] -3.97 1.42 -6.56 -2.33 5953.48 1
alpha[14] -3.56 1.18 -5.63 -2.13 6127.71 1
alpha[15] -2.85 0.94 -4.45 -1.54 7014.93 1
alpha[16] -3.64 1.19 -5.81 -2.17 5623.89 1
alpha[17] -3.67 1.24 -5.92 -2.19 5799.25 1
alpha[18] -3.91 1.35 -6.35 -2.32 5931.49 1
alpha[19] -3.38 1.04 -5.27 -2.05 5806.41 1
alpha[20] -4.01 1.47 -6.68 -2.32 5875.11 1
alpha[21] -3.30 0.99 -5.10 -2.01 5676.94 1
alpha[22] -4.11 1.54 -6.90 -2.37 6004.25 1
alpha[23] -2.96 1.00 -4.63 -1.59 6330.49 1
alpha[24] -3.83 1.27 -6.14 -2.29 5665.95 1
alpha[25] -3.94 1.37 -6.45 -2.31 5914.06 1
alpha[26] -3.37 1.03 -5.23 -2.04 5602.73 1
alpha[27] -3.67 1.19 -5.83 -2.22 5708.30 1
alpha[28] -3.22 0.98 -4.95 -1.94 5714.84 1
alpha[29] -3.76 1.26 -5.98 -2.24 5581.43 1
alpha[30] -3.65 1.18 -5.80 -2.20 6258.17 1
alpha[31] -3.30 1.02 -5.09 -1.99 6017.02 1
alpha[32] -3.76 1.23 -6.03 -2.25 6243.33 1
alpha[33] -3.34 1.04 -5.20 -2.03 5576.61 1
alpha[34] -3.62 1.18 -5.78 -2.19 5746.21 1
alpha[35] -3.68 1.17 -5.82 -2.23 5623.23 1
alpha[36] -3.15 0.97 -4.86 -1.85 6079.04 1
alpha[37] -3.86 1.32 -6.17 -2.30 5952.91 1
alpha[38] -3.87 1.35 -6.35 -2.30 5488.62 1
alpha[39] -3.56 1.13 -5.57 -2.15 5793.66 1
alpha[40] -3.60 1.17 -5.71 -2.17 6041.37 1
alpha[41] -3.69 1.16 -5.82 -2.25 5367.51 1
alpha[42] -3.63 1.15 -5.71 -2.21 5612.45 1
alpha[43] -3.73 1.26 -6.04 -2.21 5481.15 1
alpha[44] -3.48 1.10 -5.44 -2.12 5580.45 1
alpha[45] -3.43 1.09 -5.37 -2.05 5500.54 1
alpha[46] -3.04 0.97 -4.72 -1.73 6104.24 1
alpha[47] -3.24 1.00 -5.01 -1.95 5578.63 1
alpha[48] -3.96 1.40 -6.53 -2.33 5665.85 1
alpha[49] -4.11 1.55 -6.97 -2.35 6357.57 1
alpha[50] -3.19 1.05 -5.02 -1.82 6212.63 1
alpha[51] -3.57 1.12 -5.57 -2.19 5247.85 1
alpha[52] -4.02 1.42 -6.60 -2.36 5990.19 1
alpha[53] -3.39 1.09 -5.31 -2.02 5960.52 1
alpha[54] -3.30 1.00 -5.06 -2.00 5612.26 1
alpha[55] -3.66 1.25 -5.93 -2.16 5775.58 1
alpha[56] -3.84 1.33 -6.23 -2.28 5862.23 1
alpha[57] -3.46 1.17 -5.53 -2.03 6330.31 1
alpha[58] -3.20 1.05 -5.03 -1.83 5915.74 1
alpha[59] -3.43 1.04 -5.31 -2.10 5349.22 1
alpha[60] -3.58 1.12 -5.66 -2.18 5420.59 1
alpha[61] -3.58 1.20 -5.75 -2.15 6401.49 1
alpha[62] -4.27 1.67 -7.29 -2.41 6393.74 1
alpha[63] -3.58 1.16 -5.68 -2.16 5658.68 1
alpha[64] -3.31 1.04 -5.15 -1.97 5730.65 1
alpha[65] -3.48 1.06 -5.39 -2.13 5396.77 1
alpha[66] -3.42 1.04 -5.30 -2.11 5389.20 1
alpha[67] -3.62 1.16 -5.76 -2.20 5774.83 1
alpha[68] -3.44 1.04 -5.30 -2.12 5340.13 1
alpha[69] -3.55 1.13 -5.57 -2.15 5514.86 1
alpha[70] -3.39 1.09 -5.34 -2.00 5797.39 1
alpha[71] -3.55 1.20 -5.68 -2.08 5855.31 1
alpha[72] -4.25 1.65 -7.23 -2.40 5849.80 1
alpha[73] -3.24 1.07 -5.10 -1.85 6235.31 1
alpha[74] -3.60 1.16 -5.70 -2.16 5322.28 1
alpha[75] -3.46 1.08 -5.39 -2.13 5529.42 1
alpha[76] -3.64 1.17 -5.75 -2.20 5654.34 1
alpha[77] -3.25 1.04 -5.07 -1.91 5745.02 1
alpha[78] -3.57 1.17 -5.70 -2.17 5817.87 1
alpha[79] -3.85 1.30 -6.24 -2.30 5706.79 1
alpha[80] -3.32 1.07 -5.24 -1.95 5798.73 1
alpha[81] -2.75 0.97 -4.31 -1.38 7040.67 1
alpha[82] -3.03 0.91 -4.62 -1.83 6354.05 1
alpha[83] -4.15 1.52 -6.97 -2.38 5982.64 1
alpha[84] -3.21 0.95 -4.92 -1.95 5733.46 1
alpha[85] -3.81 1.28 -6.13 -2.27 5489.50 1
alpha[86] -3.80 1.28 -6.13 -2.27 5860.44 1
alpha[87] -3.79 1.31 -6.15 -2.23 6104.29 1
alpha[88] -3.37 1.07 -5.30 -1.99 5980.44 1
alpha[89] -4.23 1.64 -7.24 -2.39 6308.59 1
alpha[90] -3.49 1.08 -5.39 -2.15 5854.74 1
alpha[91] -3.73 1.25 -6.03 -2.22 5788.95 1
alpha[92] -3.60 1.17 -5.74 -2.16 5395.32 1
alpha[93] -3.91 1.38 -6.40 -2.31 5698.99 1
alpha[94] -4.19 1.56 -7.03 -2.39 5822.36 1
alpha[95] -3.83 1.32 -6.22 -2.29 5954.71 1
alpha[96] -3.38 1.01 -5.21 -2.08 5620.48 1
alpha[97] -2.85 0.91 -4.38 -1.59 6602.03 1
alpha[98] -3.42 1.08 -5.35 -2.05 5437.69 1
alpha[99] -2.73 1.03 -4.38 -1.25 7403.53 1
alpha[100] -3.83 1.30 -6.18 -2.28 5552.18 1Las transiciones divergentes persisten, pero ahora tenemos suficientes estimados de los parámetros estadísticamente independientes.
Si comparamos los modelos también podemos verificar que el mod2.3 (reparametrizado) tiene un mayor poder predictivo:
WAIC SE dWAIC dSE pWAIC weight
mod2.3 362.0922 6.547198 0.000000 NA 32.97371 0.88787193
mod2.1 367.6018 16.679799 5.509571 12.83137 23.84055 0.05648879
mod2.2 367.6321 16.242133 5.539876 12.44114 23.64734 0.05563928Veamos ahora el ajuste a los datos del modelo reparametrizado (i.e. distribución predictiva posterior). Además, verifiquemos que podemos recobrar la distribución de probabilidad generadora de los datos (\(Beta(3, 5)~+~\sigma,~\sigma \sim Normal(0, \phi)\))
post <- extract.samples(mod2.3)
est_p <- sapply(1:100, simplify = 'array', FUN =
function(i) {
# equivalente a usar link(mod2.3, dat)
i_alpha <- dat$nido[i]
i_psi <- dat$parche[i]
p <-
post$mu + post$z[, i_alpha, drop = T] + post$sigma +
post$x[, i_psi, drop = T] + post$tau +
(post$beta*dat$prop_veg[i])
inv_logit(p)
})
est <- sapply(1:100, simplify = 'array', FUN =
function(x) {
i_alpha <- dat$nido[x]
i_psi <- dat$parche[x]
p <-
post$mu + post$z[, i_alpha, drop = T] + post$sigma +
post$x[, i_psi, drop = T] + post$tau +
(post$beta*dat$prop_veg[x])
rbinom(1e3, size = dat$n[x], prob = inv_logit(p))
# equivalente a sim(mod2.3, dat)
})
par(mfrow = c(1, 2))
plot(density(p + var_parche), ylim = c(0, 2.5), ylab = 'Densidad',
xlab = expression(italic(P)['supervivencia']))
for (i in 1:100) lines(density(est_p[i, ]), lwd = 0.1)
plot(density(dat$superv), col = 'red', ylab = 'Densidad',
xlab = expression(italic('Huevos sobrevivientes')), ylim = c(0, 0.18))
for (i in 1:100) lines(density(est[i, ]), lwd = 0.2)
lines(density(dat$superv), col = 'red', lwd = 2)
En efecto el modelo está adecuadamente definido.
Comparemos ahora las distribuciones posteriores por nido y parche
par(mfrow = c(1, 2))
plot(NULL, xlim = c(0, 1), ylim = c(0, 4), ylab = 'Densidad',
xlab = expression(italic(P)['nido']))
for (i in 1:100) lines(density(inv_logit(post$z[, i])), lwd = 0.3)
plot(NULL, xlim = c(0, 1), ylim = c(0, 4), ylab = 'Densidad',
xlab = expression(italic(P)['parche']))
for (i in 1:10) lines(density(inv_logit(post$x[, i])), lwd = 0.3)
La probabilidad de supervivencia de los huevos es notablemente más variable entre nidos.
Ya sabemos que la cobertura vegetal tiene un efecto positivo sobre la supervivencia de los huevos, nosotros la simulamos, pero veamos sus implicaciones para el modelo.
v <- seq(0, 1, length.out = 100)
b_veg <- sapply(1:100, simplify = 'array', FUN =
function(i) {
inv_logit(post$beta*v[i])
})
ci_veg <- apply(b_veg, 2, quantile, c(0.025, 0.975))
plot(v, apply(b_veg, 2, mean), type = 'l', xlab = 'Cobertura de vegetación (prop)',
ylab = 'Probabilidad de supervivencia')
shade(ci_veg, v, col = col.alpha('cyan4'))
Graficamos el efecto neto de la cobertura vegetación, también podríamos estratificarlo por cada nido o parche de bosque.
4 Ejemplo
Las simulaciones estocásticas nos permite plantear el modelo generador de los datos y estudiar el desempeño de un modelo estadístico, sus supuestos y alcances, para recuperar dichos parámetros iniciales. Es sin duda ilustrativo. Ahora, apliquemos lo ya expuesto a un set de datos empírico. En este caso usaremos datos de un censo de fertilidad, específicamente el uso de anticonceptivos por mujeres bengalíes (Bangladesh) en 1988. Carguemos rethinking::bangladesh y exploremos las variables:
[1] "woman" "district" "use.contraception"
[4] "living.children" "age.centered" "urban"
[7] "district_id" d contiene 6 variables:
woman: identificador de la mujer.distric: distrito al que pertenece.use.contraception: uso (1) o no (0) de anticonceptivos.living.children: número de hijos vivos.age.centered: edad centrada.urban: contexto urbano (1) o no urbano (0). En nuestro caso lo llamaremos rural.
Nuestro objetivo será evaluar cómo el distrito, el contexto y la edad, modulan la probabilidad de que una mujer bengalí, en 1988, usara un método anticonceptivo. Definamos primero el modelo en notación matemática:
\[ \begin{align} &Uso~anticonceptivos \sim Binomial(1, P) \\ &logit(P_i) = \alpha_{urbano_i} + \theta_{distrito_i} + \beta\times edad \\ &\alpha_{urbano_i} \sim Normal(\mu, \sigma) \\ &\theta_{distrito_i} \sim Normal(0, \phi) \\ &\mu \sim Normal(0, 3) \\ &\sigma \sim exponential(1) \\ &\phi \sim exponential(1)\\ &\beta \sim Normal(0.5, 0.2) \end{align} \] Detengámonos un momento en las probabilidades previas, y discutamos mi razonamiento para su elección. Los centros urbanos en general tienen un mayor acceso a establecimientos educativos y de salud. Sería razonable esperar que la probabilidad de usar anticonceptivos sea más baja en contextos rurales comparada con los urbanos. La previa \(\mu \sim Normal(0, 3)\) implica justamente esto en la escala logit y de probabilidad, veamos:
par(mfrow = c(1, 2))
curve(dnorm(x, 0, 1.5), xlim = c(-5, 5), xlab = 'Log-odds(P)', ylab = 'Density')
plot(density(inv_logit(rnorm(1e3, 0, 1.5))), xlab = 'P', ylab = 'Density',
main = '')
Esta previa permite que el algoritmo explore, en escala logit, tanto valores negativos como positivos restringidos a \(-3\) y \(+3\). En probabilidades, implica todo el rango de probabilidades entre 0 y 1.
El contexto social, embebido en el espacial, probablemente modula la elección que hacen las mujeres. Esto es, más allá de vivir en un ambiente rural o urbano, variaciones en el sistema de creencias entre las comunidades en las que viven las mujeres, nos llevarían a esperar oscilaciones en la probabilidad de uso de anticonceptivos de acuerdo al distrito. La previa que hace justamente incorpora dicha variación eso \(\theta_{distrito_i} \sim Normal(0, \phi)\).
Finalmente, con el tiempo llegan la experiencia, madurez, la introspección (PONELE). Sea por razones éticas, ambientales, filosóficas o simplemente cansancio, es razonable suponer que la probabilidad de usar métodos anticonceptivos es mayor en mujeres más adultas. Por lo tanto, empleamos \(\beta \sim Normal(0.15, 0.2)\) una previa que permite al algoritmo explorar principalmente valores positivos de \(\beta\), pero permaneciendo relativamente conservador. Veamos:
betas_i <- rnorm(100, 0.15, 0.2)
plot(NULL, xlim = c(-10, 10), ylim = c(0, 1), xlab = 'Edad', ylab = 'Probabilidad')
for (i in 1:100) curve(inv_logit(0 + betas_i[i]*x), add = T, lwd = 0.1)
curve(inv_logit(0 + mean(betas_i)*x), add = T, lwd = 2, col = 'red')
Una vez definidas las previas, agrupemos las variables en una lista y ajustemos el modelo:
dat <-
list(
anticonc = d$use.contraception,
urbano = as.integer(ifelse(d$urban == 1, 1, 2)),
distrito = d$district_id,
edad = d$age.centered
)
me2 <-
ulam(
alist(
anticonc ~ dbinom(1, p),
logit(p) <- alpha[urbano] + theta[distrito] + beta*edad,
alpha[urbano] ~ dnorm(mu, sigma),
mu ~ dnorm(0, 3),
sigma ~ dexp(1),
theta[distrito] ~ dnorm(0, phi),
phi ~ dexp(1),
beta ~ dnorm(0.15, 2)
), data = dat, chains = 3, cores = 3, iter = 2e3, warmup = 500
) mean sd 5.5% 94.5% n_eff Rhat4
alpha[1] -0.067623135 0.119730951 -0.2613550345 0.11629924 2661.306 1.0009904
alpha[2] -0.697714646 0.089480913 -0.8412967229 -0.55744124 2550.055 1.0020527
mu -0.359706451 0.613382543 -1.2577574318 0.57912877 2277.961 1.0008869
sigma 0.803505479 0.608047107 0.2414636036 1.90131966 2221.526 0.9998454
theta[1] -0.599304600 0.202704294 -0.9227228921 -0.28170921 4759.385 1.0000798
theta[2] 0.038948422 0.335329597 -0.4934286630 0.57333786 6139.359 0.9994649
theta[3] 0.208789878 0.448681958 -0.4856305371 0.93093541 6991.348 0.9994889
theta[4] 0.277364360 0.294608936 -0.1805677048 0.74796706 6118.582 1.0001207
theta[5] 0.049472079 0.279586834 -0.3914348255 0.49617175 6168.721 1.0000162
theta[6] -0.202040920 0.247100924 -0.5984151302 0.18328980 5389.558 1.0001190
theta[7] -0.109407304 0.350072375 -0.6618394789 0.44126899 5992.130 0.9994665
theta[8] 0.101527877 0.278718386 -0.3490672465 0.53620462 5411.848 0.9997172
theta[9] -0.122873657 0.318932100 -0.6318950800 0.36903199 5737.820 1.0006100
theta[10] -0.487790376 0.388548293 -1.1130789974 0.10880375 4767.821 0.9995097
theta[11] -0.817275130 0.395063024 -1.4720090536 -0.22291037 4228.750 1.0002139
theta[12] -0.079482865 0.295096835 -0.5480543186 0.38909699 5351.155 0.9998210
theta[13] 0.064206887 0.308955223 -0.4247021249 0.55490469 5970.677 1.0004185
theta[14] 0.590036975 0.205239299 0.2725557328 0.91680047 3344.421 1.0017794
theta[15] -0.059160818 0.327432017 -0.5837396779 0.44807893 6832.870 1.0001296
theta[16] 0.434763838 0.332175877 -0.0859473266 0.97562638 5013.396 0.9995832
theta[17] -0.105885294 0.322842657 -0.6271451180 0.40292120 5891.433 0.9995406
theta[18] -0.110535507 0.260606592 -0.5306744138 0.30278628 6640.527 0.9996438
theta[19] 0.060696961 0.314500139 -0.4480547505 0.56235886 5541.327 0.9997657
theta[20] 0.124729689 0.348770116 -0.4379687217 0.68760516 6455.818 1.0000427
theta[21] -0.028042100 0.339341074 -0.5702944514 0.51639733 6082.475 0.9995438
theta[22] -0.307562262 0.357366234 -0.8998129191 0.24260176 6133.052 1.0005822
theta[23] -0.139383688 0.348576571 -0.6888358205 0.41102435 5733.781 1.0003020
theta[24] -0.496562049 0.388864112 -1.1415531687 0.09199613 4865.896 1.0004297
theta[25] 0.245790364 0.233544889 -0.1347417403 0.62060855 4715.160 1.0012794
theta[26] 0.075357514 0.360224453 -0.5067574351 0.63742546 8257.000 0.9995852
theta[27] -0.560230934 0.290611835 -1.0400880302 -0.11424822 4916.741 0.9997914
theta[28] -0.336482084 0.266137948 -0.7685811438 0.08275191 5891.932 0.9996635
theta[29] -0.227290483 0.303590220 -0.7147637386 0.25212015 6778.844 1.0005527
theta[30] 0.371483562 0.238639727 -0.0017645806 0.76108503 5257.428 1.0004287
theta[31] 0.237392829 0.284057377 -0.2159486043 0.68780935 5366.209 0.9998432
theta[32] -0.313830476 0.335727404 -0.8518602138 0.20946569 6851.961 0.9998803
theta[33] 0.021843529 0.362176793 -0.5565052465 0.59559378 8691.625 0.9993968
theta[34] 0.737474145 0.292560150 0.2854373509 1.20917261 4181.755 1.0004515
theta[35] 0.306082015 0.253072255 -0.0942720924 0.72393981 5308.822 0.9995749
theta[36] -0.012541532 0.338975127 -0.5598435193 0.52339028 6705.507 1.0000943
theta[37] 0.330857882 0.369935487 -0.2526692248 0.94058360 6546.166 1.0001492
theta[38] -0.241237493 0.367007245 -0.8302560534 0.34408371 7824.114 0.9994327
theta[39] 0.366650207 0.303106485 -0.1099512774 0.84814130 5526.012 1.0000013
theta[40] 0.050250960 0.274814925 -0.3935201019 0.49081354 5499.953 0.9998252
theta[41] 0.342069131 0.304011539 -0.1370999800 0.83395294 4520.083 0.9994641
theta[42] 0.205757053 0.369363930 -0.3890289729 0.79954716 6088.126 0.9998582
theta[43] 0.416881887 0.262705458 0.0018121456 0.84242710 5008.772 1.0003398
theta[44] -0.306242536 0.322997449 -0.8382404419 0.20666718 6958.646 0.9998799
theta[45] -0.044172478 0.282830321 -0.5013349105 0.40905733 4620.261 1.0003396
theta[46] 0.568056189 0.212349850 0.2376643239 0.91771552 3634.628 1.0003870
theta[47] 0.114518102 0.354395451 -0.4490253248 0.68650756 6586.014 0.9996146
theta[48] 0.368589873 0.265447454 -0.0519473665 0.80093663 5207.538 1.0000257
theta[49] -0.238068481 0.431643942 -0.9597736091 0.42348708 4948.379 1.0004849
theta[50] 0.224592277 0.327903030 -0.2958001919 0.75645940 5441.703 1.0000503
theta[51] 0.169842603 0.279636630 -0.2783582146 0.61936576 5002.899 1.0003238
theta[52] 0.202065160 0.234010820 -0.1580262889 0.57187789 4721.604 1.0006305
theta[53] -0.136673770 0.328904815 -0.6728171056 0.37650552 6463.956 1.0000512
theta[54] -0.318047091 0.410043974 -1.0061710008 0.30293107 5780.485 1.0001139
theta[55] 0.507416585 0.274841483 0.0744929609 0.94631989 5983.181 0.9995570
theta[56] -0.450798699 0.320002581 -0.9695122854 0.05122441 5133.193 0.9999173
theta[57] 0.154713565 0.283769109 -0.2958347699 0.61161928 6506.956 0.9998271
theta[58] -0.340627928 0.393266290 -0.9916242373 0.26366457 5405.205 0.9994447
theta[59] -0.441574047 0.323133425 -0.9653461483 0.05124577 5966.759 0.9996401
theta[60] -0.488935939 0.293576209 -0.9780675700 -0.03815926 5803.499 0.9999078
phi 0.459392769 0.080073474 0.3414493807 0.59177311 1394.408 1.0011077
beta 0.008987782 0.005349845 0.0004137012 0.01748027 7341.782 0.9995736Todos los parámetros tienen más de mil estimados estadísticamente independientes y Rhat ~ 1. Podríamos —más bien deberíamos, graficar las cadenas para corroborar que convergieron a la misma distribución posterior y exploraron el mismo espacio del parámetro. Pero, en este caso, confiaremos en la salida del modelo.
Verifiquemos el ajuste del modelo a los datos:
post2 <- extract.samples(me2)
ppcheck <- sapply(1:length(dat$anticonc), simplify = 'array', FUN =
function(x) {
urbano <- dat$urbano[x]
distrito <- dat$distrito[x]
edad <- dat$edad[x]
p <-
post2$alpha[, urbano] +
post2$theta[, distrito] +
post2$beta*edad
rbinom(1e3, 1, inv_logit(p))
})
plot(NULL, xlim = c(-0.3, 1.3), ylim = c(0, 2.8),
xlab = 'Uso de anticonceptivo', ylab = 'Densidad')
for (i in 1:200) lines(density(ppcheck[i, ]), lwd = 0.1)
lines(density(dat$anticonc), col = 'red', lwd = 2)
El modelo es confiable. Grafiquemos ahora los efectos condicionales, y verifiquemos las implicaciones del distrito y el contexto espacial en la probabilidad de que las mujeres usen algún método anticonceptivo. Primero realicemos cálculos y un poco de carpintería:
bar_urb <- unique(tibble(dist = dat$distrito,
urb = dat$urbano))
post_urb_dist <- sapply(1:nrow(bar_urb), simplify = 'array', FUN =
function(x) {
post2$alpha[, bar_urb$urb[x], drop = T] +
post2$theta[, bar_urb$dist[x], drop = T]
})
bar_urb$level <- ifelse(bar_urb$urb == 1, 'Urbano', 'Rural')
post_urb_dist <- as_tibble(post_urb_dist)
colnames(post_urb_dist) <- paste(bar_urb$level, bar_urb$dist, sep = '')
post_urb_dist <- gather(post_urb_dist)
post_urb_dist <-
post_urb_dist |>
group_by(key) |>
transmute(li = quantile(value, 0.025),
ls = quantile(value, 0.975),
mu = mean(value)) |>
ungroup() |>
unique()
post_urb_dist$dist <- gsub('^([aA-zZ]*)([0-9]{1,})', '\\2', post_urb_dist$key)
post_urb_dist$factor <- gsub('^([aA-zZ]*)([0-9]{1,})', '\\1', post_urb_dist$key)
post_urb_dist[] <- lapply(post_urb_dist, function(x) if (is.double(x)) inv_logit(x) else(x))
alpha_urb <- as_tibble(post2$alpha)
colnames(alpha_urb) <- c('Urbano', 'Rural')
contraste <- tibble(x = alpha_urb$Urbano - alpha_urb$Rural)
alpha_urb <- gather(alpha_urb)Ahora el gráfico
plot_contraste <-
ggplot(contraste, aes(x)) +
geom_density(fill = 'seagreen', color = 'seagreen', alpha = 0.5) +
labs(x = quote(italic(P['(Diferencia urbano-rural)'])), y = 'Densidad') +
geom_vline(xintercept = mean(contraste$x), linetype = 3) +
theme_bw() +
theme(panel.grid = element_blank(),
text = element_text(family = 'Times New Roman'))
plot_cond <-
post_urb_dist |>
ggplot(aes(xmin = li, xmax = ls, x = mu,
y = fct_reorder(dist, mu, .desc = T), color = factor)) +
geom_point(position = position_dodge(width = 0.3)) +
geom_linerange(position = position_dodge(width = 0.3)) +
labs(x = quote(italic(P['(usar anticonceptivos)'])), y = 'Distrito') +
theme_bw() +
theme(panel.grid = element_blank(),
legend.title = element_blank(),
legend.position = 'top',
text = element_text(family = 'Times New Roman'))
plot_contraste +
plot_cond +
plot_layout(ncol = 2)
Concluimos que las mujeres en zonas urbanas tienen, en general, un 60% más de probabilidad de usar anticonceptivos comparadas con mujeres en ambientes rurales (figura izquierda). Pero también observamos que, además del contexto espacial, dicha probabilidad está fuertemente influenciada por el distrito (figura derecha).
Ahora veamos cómo se relaciona esta probabilidad con la edad de las mujeres:
v_age <- seq(min(dat$edad), max(dat$edad), length.out = 100)
sim_edad <- sapply(v_age, simplify = 'array', FUN =
function(x) {
inv_logit(post2$beta*x)
})
sim_edad <-
as_tibble(
do.call('rbind',
apply(sim_edad, 2, simplify = 'list', FUN =
function(x) {
tibble(li = quantile(x, 0.025),
ls = quantile(x, 0.975),
mu = mean(x))
}))
)
sim_edad$x <- v_age
sim_edad |>
ggplot(aes(x = x, y = mu, ymin = li, ymax = ls)) +
geom_ribbon(alpha = 0.5, fill = 'cyan4') +
geom_line(color = 'tan1', linewidth = 1.5) +
labs(x = 'Edad (std)', y = quote(P['(usar anticonceptivos)'])) +
theme_bw() +
theme(panel.grid = element_blank(),
text = element_text(family = 'Times New Roman'))
En efecto, la edad incrementa la probabilidad de que las mujeres usen un método anticonceptivo.
Podría dedicar una sección a discutir las implicaciones de los resultados para una intervención, distrito específica, para fomentar el uso de anticonceptivos en mujeres bengalíes allá en 1988, pero no es el objetivo del ejercicio. Eso se lo dejo a quién sea que haya llegado hasta acá.
R version 4.2.1 (2022-06-23)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur ... 10.16
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] rstan_2.21.8 StanHeaders_2.21.0-7 patchwork_1.1.2
[4] magrittr_2.0.3 lubridate_1.9.2 forcats_1.0.0
[7] stringr_1.5.0 dplyr_1.1.1 purrr_1.0.1
[10] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1
[13] ggplot2_3.4.1 tidyverse_2.0.0 rethinking_2.13.2
loaded via a namespace (and not attached):
[1] Rcpp_1.0.10 mvtnorm_1.1-3 lattice_0.20-45 prettyunits_1.1.1
[5] ps_1.7.2 digest_0.6.29 utf8_1.2.3 R6_2.5.1
[9] stats4_4.2.1 evaluate_0.17 coda_0.19-4 pillar_1.8.1
[13] rlang_1.1.0 rstudioapi_0.14 callr_3.7.3 Matrix_1.5-1
[17] rmarkdown_2.17 splines_4.2.1 labeling_0.4.2 htmlwidgets_1.5.4
[21] loo_2.5.1 munsell_0.5.0 compiler_4.2.1 xfun_0.33
[25] pkgconfig_2.0.3 pkgbuild_1.4.0 shape_1.4.6 mgcv_1.8-40
[29] htmltools_0.5.3 tidyselect_1.2.0 gridExtra_2.3 codetools_0.2-18
[33] matrixStats_0.63.0 fansi_1.0.4 crayon_1.5.2 tzdb_0.3.0
[37] withr_2.5.0 MASS_7.3-57 grid_4.2.1 nlme_3.1-157
[41] DBI_1.1.3 jsonlite_1.8.4 gtable_0.3.1 lifecycle_1.0.3
[45] scales_1.2.1 RcppParallel_5.1.6 cli_3.6.0 stringi_1.7.8
[49] farver_2.1.1 ellipsis_0.3.2 generics_0.1.3 vctrs_0.6.1
[53] tools_4.2.1 glue_1.6.2 hms_1.1.2 processx_3.8.0
[57] fastmap_1.1.0 yaml_2.3.5 timechange_0.2.0 inline_0.3.19
[61] colorspace_2.1-0 knitr_1.40